#include <stdio.h>
#include <time.h>

#define abs(x) (x >= 0 ? x : -x)
#define max(x, y) (x > y ? x : y)

#define BLOCK_SIZE 5

typedef struct {
    int order;
    float * elements;
} Matrix;

__global__ void MVKernel_gm(Matrix A, Matrix X, Matrix Y0, Matrix Y1,
        int ite_max, float error, int * iter)
{
    float calc_error = 1.0;
    float value;

    int bx = blockIdx.x;
    int tx = threadIdx.x;
    // Calculate the row index of the Pd element and M
    int row = bx * BLOCK_SIZE + tx;

    if(row >= A.order) return;

    while(*iter < ite_max && calc_error > error)
    {
        (*iter)++;
        value = X.elements[row];

        for (int k = 0; k < A.order; k++)
                value -= A.elements[row * A.order + k] * Y0.elements[k];
        __syncthreads();

            Y1.elements[row] = value;
        __syncthreads();

        for(int i = 0; i < A.order; i++)
            calc_error = max(calc_error,
                (abs(Y1.elements[i] - Y0.elements[i]) / abs(Y0.elements[i])));

        for(int i = 0; i < A.order; i++)
            Y0.elements[i] = Y1.elements[i];
        __syncthreads();
    }
}

int MatMul(const Matrix A, const Matrix B, Matrix x0, Matrix x1,
	int ite_max, float error)
{
    int iter = 0;
    int * d_iter;
    cudaMalloc(&d_iter, sizeof(int));
    cudaMemcpy(d_iter, &iter, sizeof(int), cudaMemcpyHostToDevice);

    Matrix d_A;
    d_A.order = A.order;
    size_t size = A.order * A.order * sizeof(float);
    cudaMalloc(&d_A.elements, size);
    cudaMemcpy(d_A.elements, A.elements, size, cudaMemcpyHostToDevice);

    Matrix d_B;
    d_B.order = B.order;
    size = B.order * A.order * sizeof(float);
    cudaMalloc(&d_B.elements, size);
    cudaMemcpy(d_B.elements, B.elements, size, cudaMemcpyHostToDevice);

    Matrix d_x0;
    d_x0.order = x0.order;
    size = x0.order * A.order * sizeof(float);
    cudaMalloc(&d_x0.elements, size);
    cudaMemcpy(d_x0.elements, x0.elements, size, cudaMemcpyHostToDevice);

    Matrix d_x1;
    d_x1.order = x1.order;
    size = x1.order * A.order * sizeof(float);
    cudaMalloc(&d_x1.elements, size);

    dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);
    dim3 dimGrid((A.order + dimBlock.x - 1) / dimBlock.x,
            (A.order + dimBlock.y - 1) / dimBlock.y);

    /*for(int i = 0; i < A.order; i++)
        printf("%f ", x0.elements[i]);
    printf("\n");*/

    MVKernel_gm<<<dimGrid, dimBlock>>>(d_A, d_B, d_x0, d_x1,
            ite_max, error, d_iter);
    cudaThreadSynchronize();

    cudaMemcpy(x1.elements, d_x1.elements, size, cudaMemcpyDeviceToHost);
    cudaMemcpy(&iter, d_iter, sizeof(int), cudaMemcpyDeviceToHost);

    /*for(int i = 0; i < A.order; i++)
        printf("%f ", x1.elements[i]);
    printf("\n");*/

    cudaFree(&d_iter);
    cudaFree(d_A.elements);
    cudaFree(d_B.elements);
    cudaFree(d_x0.elements);
    cudaFree(d_x1.elements);

    return iter;
}


int main(int argc, char * argv[])
{
    FILE *ifp; /* Variavel que contem o ponteiro para arquivo de entrada */

    int order = 0; /* Variavel que indica a ordem da matriz */
    int row_test = 0; /* Variavel que indica a linha a ser avaliada */
    int i; /* Variavel contadora */
    int j; /* Variavel contadora */
    int ite_max = 0; /* Variavel que indica o numero maximo de iteracoes */
    int iter = 0; /* Variavel que indica o numero de iteracoes */

    float error = 0.0; /* Variavel que indica erro para a condicao de parada */
    float res = 0; /* Variavel que contem o resultado final */
    float maxt = 0.0;
    float maxl = 0.0;

    /* Matrizes da CPU ou host */
    Matrix x0; /* Valores dos x0 */
    Matrix x1; /* Valores dos x1 */
    Matrix mat_a; /* Variavel que contem a matriz A */
    Matrix mat_a_e; /* Variavel que contem a matriz A estrela */
    Matrix mat_b; /* Variavel que contem a matriz B */
    Matrix mat_b_e; /* Variavel que contem a matriz B estrela */

    /* Verifica o uso correto do programa */
    if(argc != 2)
    {
        printf("Usage: ./sequencial [inputfile]\n");
        return 1;
    }

    /* Verifica se foi possivel abrir o arquivo */
    if((ifp = fopen(argv[1], "r")) == NULL)
    {
        printf("Could not open the file %s\n", argv[1]);
    }

    /* Recebe valores de definicao da matriz */
    fscanf(ifp, "%d\n%d\n%f\n%d", &order, &row_test, &error, &ite_max);

    /************************************************/
    /* Alocacao de memoria para a CPU ***************/
    /************************************************/
    /* Aloca memoria para a matriz A */
    mat_a.elements = (float *) malloc(order * order * sizeof(float));
    /* Aloca memoria para a matriz A estrela */
    mat_a_e.elements = (float *) malloc(order * order * sizeof(float));
    /* Aloca memoria para a matriz B */
    mat_b.elements = (float *) malloc(order * sizeof(float));
    /* Aloca memoria para a matriz B estrela */
    mat_b_e.elements = (float *) malloc(order * sizeof(float));
    /* Aloca memoria para os valores de x0 */
    x0.elements = (float *) malloc(order * sizeof(float));
    /* Aloca memoria para os valores de x1 */
    x1.elements = (float *) malloc(order * sizeof(float));

    /* Recebe os valores do arquivo para a matriz A */
    for(i = 0; i < order; i++)
        for(j = 0; j < order; j++)
            fscanf(ifp, "%f", &mat_a.elements[(i * order) + j]);

    /* Recebe os valores do arquivo para a matriz B */
    for(i = 0; i < order; i++)
        fscanf(ifp, "%f", &mat_b.elements[i]);

    /* Calcula o criterio de convergencia e normaliza a matriz A e B */
    for(i = 0; i < order; i++)
    {
        maxl = 0;
        for(j = 0; j < order; j++)
        {
            if(i != j)
            {
                mat_a_e.elements[(i * order) + j] =
                    mat_a.elements[(i * order) + j] /
                    mat_a.elements[(i * order) + i];
                maxl += mat_a_e.elements[(i * order) + j];
            }
            else mat_a_e.elements[(i * order) + j] = 0;
            maxt = max(maxt, maxl);
        }

        mat_b_e.elements[i] = mat_b.elements[i] /
            mat_a.elements[(i * order) + i];
        x0.elements[i] = mat_b_e.elements[i];
    }

    if(maxt > 1)
        printf("Convergence criterion not satisfied.\n");
    else
        printf("Convergence criterion satisfied.\n");

    mat_a_e.order = order;
    mat_b_e.order = 1;
    x0.order = 1;
    x1.order = 1;

    iter = MatMul(mat_a_e, mat_b_e, x0, x1, ite_max, error);

    /* Efetua o calculo do resultado */
    for(i = 0; i < order; i++)
        res += mat_a.elements[(row_test * order) + i] * x1.elements[i];

    /* Impressao da saida */
    printf("-----------------------------------------\n");
    printf("Iterations: %d\n", iter);
    printf("RowTest: %d => [%f] =? [%f]\n",
            row_test, res, mat_b.elements[row_test]);
    printf("-----------------------------------------\n");

    /* Limpa a memoria alocada no host e o arquivo aberto */
    free(mat_a.elements);
    free(mat_a_e.elements);
    free(mat_b.elements);
    free(mat_b_e.elements);
    free(x0.elements);
    free(x1.elements);
    fclose(ifp);
    return 0;
}
